module func_int_test; ! integration
use numeric_quadr;
use numeric_deriv
use func_factory;
implicit none;

contains 

logical function func_int_test_all() result(ok);  implicit none; 
  call dbg_info(">"); 
  call dbg_info("->func_int_test_all()"); 
  ok = .TRUE.
  ok = new_func_int_test();
  return
end function

logical function new_func_int_test() result(ok);  implicit none; 
  integer :: nv;
  real(R14) :: minv, maxv, res, res2, MAX_DRV_ERR = 1.d-18;
  real(R14), dimension(:), pointer :: px, py, py2, pdy, pf, xq;      
  call dbg_info("-->new_func_int_test()"); 
  ok = .TRUE.

  nv = 5;  minv = D0; maxv = D1;
  !px => new_step_grid(minv, maxv, nv);
  px => new_vec(nv);  
  call load_step_grid(px, 1, nv, minv, maxv);
  call assert_max_zero("  5.px(1)=minv", px(1) - minv);
  call assert_max_zero("  5.px(2)=minv", px(2) - (maxv+minv)/4.);
  call assert_max_zero("  5.px(nv)=maxv", px(nv) - maxv);
  
  xq => new_step_qdr(px); ! x-quadratures
  call dbg_vec("xq=", xq);
  call assert_max_zero("  xq(1)", xq(1) - D05*7D0  /45D0);
  call assert_max_zero("  xq(2)", xq(2) - D05*32D0 /45D0);
  call assert_max_zero("  xq(3)", xq(3) - D05*12D0 /45D0);
  call assert_max_zero("  xq(4)", xq(4) - D05*32D0 /45D0);
  call assert_max_zero("  xq(5)", xq(5) - D05*7D0  /45D0);

  py => new_vec_polynom(px, (/D1, D0, D0/)); 
  res = calc_int(xq, py); 
  call assert_max_diff("int 1.=", res, 1D0, ERR15);

  ! new x-grid-9
  call del_vec3(px, py, xq);
!  px => new_step_grid(D0, D1, 9);
  px => new_vec(9);  
  call load_step_grid(px, 1, 9, D0, D1);
  py => new_vec_polynom(px, (/D1, D0, D0/)); 
  xq => new_step_qdr(px); ! x-quadratures
  res = calc_int(xq, py); 
  call assert_max_diff("9.int 1.=", res, D1, ERR15);
  
  ! new x-grid-13
  call del_vec3(px, py, xq);
!  px => new_step_grid(D0, D1, 13);
  px => new_vec(13);  
  call load_step_grid(px, 1, 13, D0, D1);
  py => new_vec_polynom(px, (/D1, D0, D0/)); 
  xq => new_step_qdr(px); ! x-quadratures
  res = calc_int(xq, py); 
  call assert_max_diff("13.int 1.=", res, D1, ERR15);
  
  call del_vec(py);
  py => new_vec_polynom(px, (/D0, D1/)); 
  res = calc_int(xq, py); 
  res2 = D05;
  call assert_max_diff("13.x^2=", res, res2, ERR15);

  call del_vec(py);
  py => new_vec_polynom(px, (/D0, D0, D1/)); 
  res = calc_int(xq, py); 
  res2 = D1/D3;
  call assert_max_diff("13.x^3=", res, res2, ERR15);

  call del_vec(py);
  py => new_vec_polynom(px, (/D1, -D1, D1, -D1, D1/)); 
  res = calc_int(xq, py); 
  res2 = D1 - D05 + D1/D3 - D1/D4 + D1/D5;
  call assert_max_diff("13.1-1+1-1+1=", res, res2, ERR15);

  ! new x-grid-21
  call del_vec3(px, py, xq);
  !px => new_step_grid(D0, PI, 21);
  px => new_vec(21);  
  call load_step_grid(px, 1, 21, D0, PI);

  xq => new_step_qdr(px); ! x-quadratures
  py => new_func_sin(px);
  res = calc_int(xq, py); 
  res2 = D2;
  call assert_max_diff("21.sinx=", res, res2, ERR7);

  py = COS(px);
  res = calc_int(xq, py); 
  res2 = D0;
  call assert_max_diff("21.(1-cosx)=", res, res2, ERR15);

  ! new x-grid-301
  call del_vec3(px, py, xq);
!  px => new_step_grid(D0, PI, 301);
  px => new_vec(301);  
  call load_step_grid(px, 1, 301, D0, PI);
  xq => new_step_qdr(px); ! x-quadratures
  py => new_func_sin(px);
  res = calc_int(xq, py); 
  res2 = D2;
  call assert_max_diff("301.sinx=", res, res2, ERR14);

  py = D1-COS(px);
  res = calc_int(xq, py); 
  res2 = PI;
  call assert_max_diff("21.(1-cosx)=", res, res2, ERR14);

  ! new x-grid-101
  call del_vec3(px, py, xq);
!  px => new_step_grid(-D1, D1, 101);
  px => new_vec(101);  
  call load_step_grid(px, 1, 101, -D1, D1);
  xq => new_step_qdr(px); ! x-quadratures
  py => new_func_exp(px);
  res = calc_int(xq, py); 
  res2 = EXP(D1)-EXP(-D1);
  call assert_max_diff("101.exp=", res, res2, D4*ERR13);

  call del_vec3(px, py, xq);
  return
end function

end module 